Semiclassical theory for spatial density oscillations in fermionic systems 
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We investigate the particle and kinetic-energy densities for a system of A'^ fermions bound in a 
local (mean-field) potential V{r). We generalize a recently developed semiclassical theory [J. Roccia 
and M. Brack, Phys. Rev. Lett. 100, 200408 (2008)], in which the densities are calculated in terms 
of the closed orbits of the corresponding classical system, to D > 1 dimensions. We regularize the 
semiclassical results (i) for the U(l) symmetry breaking occurring for spherical systems at r = 
and (ii) near the classical turning points where the Friedel oscillations are predominant and well 
reproduced by the shortest orbit going from r to the closest turning point and back. For systems with 
spherical symmetry, we show that there exist two types of oscillations which can be attributed to 
radial and non-radial orbits, respectively. The semiclassical theory is tested against exact quantum- 
mechanical calculations for a variety of model potentials. We find a very good overall numerical 
agreement between semiclassical and exact numerical densities even for moderate particle numbers 
A*'. Using a "local virial theorem", shown to be valid (except for a small region around the classical 
turning points) for arbitrary local potentials, we can prove that the Thomas- Fermi functional ttp [p] 
reproduces the oscillations in the quantum-mechanical densities to first order in the oscillating parts. 

PACS numbers: 03.65.Sq, 03.75.Ss, 05.30.Fk, 71.10.-w 



I. INTRODUCTION 

Recent experimental success confining fermion gases 
in magnetic traps Q has led to renewed interest in the- 
oretical studies of confined degenerate fermion systems 
at zerojljl, i, i, i, 0, i, B M, [H and finite tempera- 
tures IJ. According to the density functional theory 
(DFT)[linil,[ill, the local particle density p(r) is the 
key ingredient of a system of interacting fermions in that 
it contains all information about its ground state. In 
this paper we study the oscillations in the particle den- 
sity yo(r) and in different forms of the kinetic-energy den- 
sity of N fermions bound in a local potential V{r). Al- 
though wc treat the particles as non-interacting, wc keep 
in mind that this potential models the self-consistent 
Kohn-Sham (KS) potential obtained for an interacting 
system in the mean-field approximation. We shall also 
consider potentials with infinitely steep walls, so-called 
"billiards" , which have been shown to be good approxi- 
mations to the self-consisten mean fields of quantum dots 
pTf or metal clusters [l^ with many particles. 

A semiclassical theory for spatial density oscillations 
has been developed recently in [l^. Using Gutzwiller's 
semiclassical Green function [2^, expressions for the os- 
cillating parts of spatial densities of fermionic systems 
were given in terms of the closed orbits of the correspond- 
ing classical system. The semiclassical theory was shown 
in to reproduce very accurately the quantum oscil- 
lations in the spatial densities of one-dimensional sys- 
tems, even for moderate particle numbers N, and some 
general results have also been given for arbitrary higher- 
dimensional spherical potentials V{r). 

In this paper, we present in more detail the semiclas- 
sical closed-orbit theory developed in [lH and apply it 
explicitly for a variety of potentials in D > 1 dimensions. 
We find overall a good agreement between the quantum- 
mechanical and the semiclassical densities. 



The paper is organized as follows. In Sec.|TT]we give the 
basic definitions of the quantum-mechanical spatial den- 
sities. In Sec. Ill BI we discuss the asymptotic (extended) 
Thomas-Fermi (TF) limits for N ^ oo and emphasize 
the existence of two types of density oscillations occur- 
ring in potentials for D > 1 with spherical symmetry 
(except for isotropic harmonic oscillators). 

Sec. lIIIl is devoted to the semiclassical closed-orbit the- 
ory for spatial density oscillations. In Sees. IIIlXl - IIIIDl 
we review the basic equations and former results, in- 
cluding also details that were not presented in [l^. In 
Sees, nil El and IIII Fl wc extend the semiclassical theory 
to higher-dimensional systems {D > 1) and test its re- 
sults for various model potentials against exact quantum- 
mechanical densities. In Sect. IIII E3] we discuss the reg- 
ularization necessary in spherical systems for D > I near 
the center (r = 0), where a U(l) symmetry breaking oc- 
curs for r > 0. In a separate publication [2l|, wc have 
presented the analytical determination and classification 
of all closed orbits in the two-dimensional circular billiard 
and give analytical results of the semiclassical theory for 
the spatial density oscillations in this system. Some of 
the numerical results for the densities are included in Sec. 
IIII E 41 of the present paper. 

In Sec. lIVI we present regularizations of the spatial den- 
sities near the classical turning points, where the semi- 
classical theory diverges, both for smooth potentials and 
for billiard systems. 

Sec. |V] contains some general results valid for finite 
fermion systems such as trapped fermionic gases or 
metallic clusters. We discuss there, in particular, a "lo- 
cal virial theorem" and, as its direct consequence, the 
extended validity of the TF functional ttf [p] ■ 

Throughout this paper, wc only treat the zero-tem- 
perature ground state of an TV-particle system. In the 
Appendix o we outline how to include finite temperatures 
for grand-canonical ensembles in the semiclassical theory. 
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II. THE QUANTUM-MECHANICAL DENSITIES 
A. Basic definitions and ingredients 

Let us recall some basic quantum-meclianical defini- 
tions, using the same notation as in [l0| . We start from 
the stationary Schrodinger equation for particles with 
mass m, bound by a local potential V{r) with a discrete 
energy spectrum {En}: 



-— V2 + y(r) [>0„(r) (r). (1) 



We order the spectrum and choose the energy scale such 

that < El < E2 < ■ ■ ■ < En < We consider a 

system with an even number TV of fermions with spin 
s = 1/2 filling the lowest levels, and define the particle 
density by 



p(r) 2^|0„(r)p, 



J p{r) d^r = iV. 



(2) 



Hereby A is the Fermi energy and the factor 2 accounts 
for the fact that due to spin and time-reversal symme- 
try, each state n is at least two-fold degenerate. Further 
degeneracies, which may arise for D > 1, will not be 
spelled out but included in the summations over n. For 
the kinetic-energy density, we consider two different def- 
initions [231 



r(r) 2 J] C(r)V2</)„(r) , (3) 



n(r) 



2m 



(4) 



£;„<A 



which upon integration both yield the exact total kinetic 
energy. Due to the assumed time-reversal symmetry, the 
two above functions are related by 



1 fi^ 

T(r) =Ti(r)--— VV(r). 

2 Zm 



(5) 



An interesting, and for the following discussion conve- 
nient quantity is their average 



e(r) 



1 



'(r)+ri(r)]. 



(6) 



For harmonic oscillators it has been observed 0, [13, H^l 
that inside the system (i.e., sufRciently far from the sur- 
face region), ^(r) is a smooth function of the coordinates, 
whereas T(r) and ri(r), like the density p{r), exhibit 
characteristic shell oscillations that are opposite in phase 
for T and ri. We can express r(r) and ri(r) in terms of 
C(r) and V2p(r): 



n(r) = e(r) + i|^VVr): 



(7) 
(8) 



so that p(r) and ^ (r) can be considered as the basic den- 
sities characterizing our systems. Eqs. ([2]) - ([8]) are exact 
for arbitrary potentials V{r). For any even number 
of particles they can be computed once the quantum- 
mechanical wave functions 0Ti(r) are known. As men- 
tioned in the introduction, the potential V{r) can be 
considered to represent the self-consistent mean field of 
an interacting system of fermions obtained in the DFT 
approach. The single-particle wavefunctions (^Ti(r) are 
then the Kohn-Sham orbitals [l^ and p{r) is (ideally) 
the ground-state particle density of the interacting sys- 
tem. 

For later reference we express the densities ([2])- ([4]) in 
terms of the Green function in the energy representation, 
which in the basis {(/)„(r)} is given by 



G{E, r,r') 



^ E 



En 



(e > 0) . (9) 



Using the identity 1/{E + ie- En) = V[l/iE - En)] - 
iTTS{E — En), where V is the Cauchy principal value, one 
can write the densities as 

p(r) = --Im dEGiE,r,r%,^r, (10) 

Jo 

T(r) = -^Im/ di?V2,G(i?,r,r')|r'=r, (11) 
27rm Jq 

n(r) = Im / di;VrVr'G(£;,r,r')|r'=r, (12) 

27rm Jn 



whereby the subscript of the nabla operator V denotes 
the variable on which it acts. 

The density of states g{E) of the system ([T]) is given by 
a sum of Dirac delta functions, which can be expressed 
as a trace integral of the Green function: 

g{E) = Y,S{E-En)^-^lm J d^r G(£;, r, r')|r'=r • 



The particle number can then also be obtained as 



N = NiX) = 2 / dEg{E) 



(13) 



(14) 



Due to the discreteness of the spectrum, N{X) is a 
monotonously increasing staircase-function and conse- 
quently the function X{N), too, is a monotonously in- 
creasing staircase-function. 



B. Asymptotic quantum- mechanical results 

1. Thomas-Fermi limits and oscillating parts 

In the limit N ^ 00, the densities are expected to go 
over into the approximations obtained in the Thomas- 
Fermi (TF) theory 0|. These are given, for any local 
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potential ^(r), by 

A 1 

PTF{r) 



(^r^ [Atf - Vir)] 



D r(f ) \2iTh?) 

(Ti)TF(r) = ^TF(r) = TTF(r), 



D/2 



(15) 
(16) 



TTF(r) 



1 / m \^/2 



{D+2) r(f ) \2^h^) 



[Atf - V{v)] 



D/2+1 



(17) 

These densities are defined only in the classically allowed 
regions where Atf > ^(r), and the Fermi energy Atf 
is defined such as to yield the correct particle number 

upon integration of pTF(r) over all space. The direct 
proof that the quantum- mechanical densities, as defined 
in Sec. |TT] in terms of the wavefunctions of a smooth po- 
tential, reach the above TF limits for ^ oo is by no 
means trivial. It has been given for isotropic harmonic 
oscillators in arbitrary dimensions in Ref . [lO| . 

The TF densities ([ini)-(II71) fulfill the following func- 
tional relation: 



TTF(r) = TTF[PTF(r)] 



AttD 



2m {D + 2) 



4 \ 2 



-|2/r> 



4+^/^(r), (18) 



which will be investigated further below. 

For smooth potentials in D > 1 dimensions, next-to- 
leading order terms in 1/N modify the smooth parts of 
the spatial densities, which are obtained in the extended 
Thomas-Fermi (ETF) model as corrections of higher or- 
der in h through an expansion in terms of gradients of 
the potential [2^. These corrections usually diverge at 
the classical turning points and can only be used suffi- 
ciently far from the turning points, i.e., in the interior 
of the system. We do not reproduce the ETF densities 
here but refer to [2^ (chapter 4) where they are given for 
arbitrary smooth potentials in _D = 2 and 3 dimensions, 
and to [lO[ where explicit results are given for spherical 
harmonic oscillators in Z? = 2 and 4 dimensions. 

This leads us to decompose the densities in the follow- 
ing way: 



P(i') = p(E)TF(r) + '^p(r), 

T(r) = r(E)TF(r) + ^T(r), 

Ti(r) = (ri)(E)TF(r) + '5-ri(r) 

^(r) = e(E)TF(r) + <5e(r). 



(19) 
(20) 
(21) 
(22) 



For D ~ 1 and for billiard systems [27[ , the subscripts TF 
and the explicit relations (jTS]) - p7|) hold. The oscillating 
parts Sp(r) etc. are the main objects of this paper. 



2. Two types of oscillating parts in spherical systems 

We have investigated the density oscillations in vari- 
ous potentials in Z? > 1 dimensions with radial symmetry 
such that V{r) — V{r), where r = |r|. We found that. 



generally, there exist two types of oscillations in their 
spatial densities: 

(i) regular, short-ranged oscillations with a constant 
wavelength in the radial variable r over the whole re- 
gion, and 

(ii) irregular, long-ranged oscillations whose wavelength 
decreases with increasing r. 

An example is shown in Fig. [T] for a spherical billiard 
with unit radius containing N ~ 100068 particles. Note 
the irregular, long-ranged oscillations of <^(r) around its 
bulk value [23| ^tf seen in the upper panel. In the lower 
panel, where we exhibit only an enlarged region around 
the bulk value, we see that T(r) and Ti(r) oscillate reg- 
ularly around £,(r), but much faster than ^(r) itself and 
with opposite phases. The same two types of oscillations 
are also found in the particle density p{r). 



0.6 



0. 

^ 0.5 
^ — ^ 

0.4 

§ 0.58 
^ 0.56 



0.54 




0.0 0.1 0.2 0.3 0.4 



0.5 
r/R 



0.6 0.7 0.8 0.9 1.0 



FIG. 1: Kinetic-energy density profiles of a 3-D spherical bil- 
liard with N = 100068 particles (units: h'^/2m = i? = 1; all 
densities are divided by A''^''^). Upper panel: ^(r) (solid line) 
and its constant TF value ^tp (dashed). Lower panel: T(r) 
(dashed), ri(r) (dotted) and ^(r) (solid line). Note that in 
both panels, the vertical scale does not start at zero. 

For radial systems, we can thus decompose the oscil- 
lating parts of the spatial densities defined in (fTO|) 
as follows: 



5p{r) = Srp{r) + 5\„p{r) , 
5r(r) = (SrT(r) -f- 5{„T{r) , 

5Ti{r) = (SrTl(r) + (5irrTl(r) , 



(23) 
(24) 
(25) 
(26) 



Here the subscript "r" denotes the regular, short-ranged 
parts of the oscillations, while their long-ranged, irregular 
parts are denoted by the subscript "irr". We emphasize 
that this separation of the oscillating parts does not hold 
close to the classical turning points. 

As we see in Fig. [T] and in later examples, the oscillat- 
ing parts defined above fulfill the following properties in 
the interior of the system (i.e., except for a small region 
around the classical turning points): 
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a) For D > 1, the irregular oscillating parts of r(r) and 
Ti(r) are asymptotically identical and equal to S^{r): 



(27) 



b) The irregular oscillations are absent (i.e., asymptoti- 
cally zero) in the densities of all potentials in D = 1 and, 
also, in isotropic harmonic oscillators (see Ref. [l^) and 
in linear potentials (see [1^) for arbitrary D. 

and Ti(r) are 



c) The regular oscillating parts of T(r) 
asymptotically equal with opposite sign: 



(5rr(r) 



-drTi{r) . 



(28) 



(This relation holds in particular for isotropic harmonic 
oscillators, for which it has been derived [lO| asymptoti- 
cally for N ^ oo from quantum mechanics.) 

These numerical findings will be understood and ex- 
plained within the semiclassical theory developed in the 
following. 

Henceforth, the symbol 6 will always denote the sum 
of both types of oscillating parts and the subscripts will 
only be used if reference is made to one particular type 
of oscillations. 



III. SEMICLASSICAL CLOSED-ORBIT 
THEORY 

In this section we present the semiclassical theory, initi- 
ated by Gutzwiller (see [2^ and earlier references quoted 
therein, and [20|). for the approximate description of 
quantum oscillations in terms of classical orbits. In Sec. 
nil Al we recall the trace formula for the density of states, 
and in Sec. IIIIBI we present the newly developed the- 
ory for spatial density oscillations In both cases, 
we limit ourselves - as in the previous section - to 
non-interacting fermions in a local potential V{r). The 
inclusion of finite temperatures in the semiclassical the- 
ory is dealt with in Appendix □ 



A. Brief review of periodic orbit theory for the 
density of states 

Before deriving semiclassical expressions for the spatial 
densities, we remind the reader of the periodic orbit the- 
ory (POT) for the density of states. The starting point 
is the semiclassical approximation of the Green function 
^ which was derived by Gutzwiller [20| : 



(29) 



The sum runs over all classical trajectories 7 leading from 
a point r to the point r' at fixed energy E. S^{E, r,r') 
is the action integral taken along the trajectory 7 



S^{E,r,r')^ p(i?,q).dq. 



(30) 



whereby p{E, r) is the classical momentum 



piE,v)^-^2m[E-Vir)] 



(31) 



T)^ is the 



defined only inside the classically allowed region where 
E > V{r); its modulus is denoted by p{E,r 
Van Vleck determinant: 

''^ = piE,r)piE,r') ''-' P. = det(ap,/ar',) , 

(32) 

where p_L and are the initial momentum and final 
coordinate, respectively, transverse to the orbit 7. The 
Morse index fi^ counts the sign changes of the eigenval- 
ues of the Van Vleck determinant along the trajectory 7 
between the points r and r'; it is equal to the number of 
conjugate points along the trajectory [33|. The prefactor 
in (|^ is given by 



ao ~ 27r(2i7rfi) 



-{D+l)/2 



(33) 



The approximation (j29p of the Green function is now 
inserted into the r.h.s of (fT3|) for the density of states 
g{E). Since r' = r in the trace integral of only 
closed orbits contribute to it. The running time Tj{E, r) 
of these orbits, i.e., the time it takes the classical particle 
to run though the closed orbit, is given by 



T^iE,r) 



dS^{E, 



dE 



(34) 



It was shown by Berry and Mount [3l| that to leading 
order in H, the orbits with zero running time, T^{E, r) = 
0, yield the smooth TF value of g{E). In systems with 
D > 1 higher-order terms in h also contribute, which can 
also be obtained from the ETF model (see, e.g., chapter 
4 of [1^). Separating smooth and oscillatory parts of the 
density of states by defining 



giE) := g{E) + Sg{E) 



(35) 



the oscillating part Sg{E) is, to leading order in h, given 
by the semiclassical trace formula 



SgiE)~J2^PoiE) cos 



PO 



1 TT 

- SpoiE) — — (7po 



(36) 



where the sum runs over all periodic orbits (POs). For 
systems in which all orbits are isolated in phase space, 
Gutzwiller [29j derived explicit expressions for the am- 
plitudes Ay>o{E), which depend on the stability of the 
orbits, and for the Maslov indices crpo- Performing the 
trace integral in (|13p along all directions transverse to 
each orbit 7 in the stationary phase approximation (SPA) 
leads immediately to the periodicity of the contributing 
orbits. The Maslov index crpo collects all phases occur- 
ring in ([^^ and in the SPA for the trace integral (see 32 1 
for detailed computations of (Tpo). It has been shown 33 1 
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that (Tpo is a canonical and topological invariant prop- 
erty of any PO. Spo{E) is the closed action integral 



SpoiE) = / p(£;,q).dq. 

Jpo 



(37) 



For smooth one-dimensional potentials, the trace formula 
is particularly simple and reads 



(D-- 



fc=i 



(38) 



where the sum is over the repetitions fc > 1 of the primi- 
tive orbit with action Si{E) and period Ti{E) = S[{E). 
Equation (|38p is equivalent to the sum of delta functions 
in (fT3|) , using the spectrum obtained in the WKB approx- 
imation ^26. . .S^ . For systems with D > 1 with continu- 
ous symmetries (and hence also for integrable systems), 
the same type of trace formula (|36p holds, but the sum- 
mation includes all degenerate families of periodic orbits 
and the amplitudes „4po (E) and indices crpo have differ- 
ent forms. For an overview of various trace formulae and 
the pertinent literature, as well as many applications of 
the POT, we refer to [H. 



Semiclassical approximation to the spatial 
densities 



In order to derive semiclassical expressions for the spa- 
tial densities defined in Sec. [Til we start from the expres- 
sions given in the equations (fTO)) - p^ . which are func- 
tions of r and the Fermi energy A, and replace the exact 
Green function G{E, r, r') by its semiclassical expansion 
(j29|) . The energy integration can be done by parts, using 
((32)) and (p4)) . and to leading order in h we obtain for the 
particle density 



p(A,r) 



2mh 



7rp(A,r)^*'"-5 T,(A,r) 



i piS'T(A,r,r)-i^i^ 



(39) 

Again, the orbits with zero running time T[E, r) = 
yield, to leading order in ft, the smooth TF particle den- 
sity ; the proof given in [3l| for the density of states 
applies also to the spatial densities discussed here. Like 
for the density of states, higher-order h corrections con- 
tribute also to the smooth part of p(r) in Z? > 1 and will 
be included in their ETF expressions. The periodic or- 
bits (POs), too, can only contribute to the smooth part 
of p(r), since their action integrals ([57)1 are independent 
of r and hence the phase in the exponent of (|39p is con- 
stant. Thus, a priori only non-periodic orbits (NPOs) 
contribute to the oscillating part of /o(r). The same holds 
also for the other spatial densities, so that we can write 



their semiclassical approximations as [IE 
R-C 2^ 



5p{v) 
5t{t) 



7rp(A,r; 
hp(X,r) 



NPO 



r(A,r) 



Re a„ 



NPO 



T(A,r) 



p»*(A,r) 



(40) 



(41) 



„»*(A,r) 



<5ri(r) ~ Re > Q(A, r) ~ \ 

TT T(A,r) 

(42) 

The sums are only over non-periodic orbits (NPOs) that 
lead from a point r back to the same point r. For con- 
venience, we have omitted the subscript "NPO" from all 
quantities in the above equations. The phase function 
<I'(A, r) is given by 



$(A,r) 



S{X,r,r)/h-fi-. 



(43) 



The quantity Q(A, r) appearing in (|42|) for (5Ti(r) is de- 
fined as 



Q(A,r) 



p(A,r).p(A,r') 



p2(A,r) 



cos[0(p,p')], (44) 



where p and p' are the short notations for the initial and 
final momentum, respectively, of a given closed orbit 7 
at the point r. These are obtained also from the action 
integral ([50)) by the canonical relations 



Vr5-,(A,r, r') 



-p, Vr'5^(A,r,r') 



, = P ■ 

(45) 

Since Q in ([^1)1 depends on the angle 6 between p and 
p', it may be called the "momentum mismatch function" , 
being -1-1 for p = p' (i.e., for POs) and —1 for p = — p' 
(e.g., for self-retracing NPOs). 

Note that the upper limit A of the energy integral in 
([To)) - ([T2)) has been replaced here by the smooth Fermi 
energy A defined by 



N = 2 I AEg{E) 
/o 



A = A + (5A . 



(46) 



The reason for this is the following. Since A(A^) is a non- 
smooth staircase function, as mentioned at the end of 
Sec. [m it is natural to expand it around its smooth part 
A which can be identified with its TF value Atf (or Aetf 
for D > 1). Taylor expanding equation ([T^ using ([55)1 
up to first order in 5\, we easily obtain an expression for 
its oscillating part (cf. [H]): 



(5A: 



1 



5'etf(A) J a 



dE6g{E) . 



(47) 



The quantity 6X is of higher order in ft than A and can be 
considered as a small semiclassical correction; the Sg{E) 
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in the integrand may be expressed through the trace for- 
mula ([55]) . Now, the contribution of the zero- length or- 
bits to ([39|) yields formally the smooth (E)TF density, 
but taken at the exact (quantum) value of A. The den- 
sity should therefore be developed around the smooth 
(E)TF value A before it can be identified with the stan- 
dard (E)TF density. Its first variation with 6X leads to a 
further smooth contribution which should be taken into 
account. The same holds for the other densities. The 
contribution of all finite-length orbits to ([5^ is of higher 
order in h than the leading smooth (ETF) terms, so it is 
consistent to evaluate them at A. 

In one-dimensional systems, all smooth terms can be 
exactly controlled. The smooth part of the density may 
be written as 



Ptf(A, x) ~ /9tf(A, x) + 6X 



dpTF(A, a;) 



dA 



(48) 



The first term on the r.h.s. is the standard TF density 
(for D = 1). The second term, using (|38p and the fact 
that gTF(ATF) = T'i(ATF)/27r?i for D = 1, is found to ex- 
actly cancel the contribution of the periodic orbits to pQ]) 
(evaluated at A), which has been explicitly calculated in 
jigj and given in Eq. (22) there. 

For D > 1 dimensions, we cannot prove that the same 
cancellation of smooth terms tak es p lace. Furthermore, 
for the circular billiard treated in [21| it is shown that the 
contributions of periodic and non-periodic orbits cannot 
be separated in the vicinity of bifurcations that occur for 
D > 1 under variation of r. For arbitrary local potentials 
in D > I dimensions, it is in general a difficult task to 
evaluate all nonperiodic closed orbits. In non-integrable 
systems, the number of POs is known to grow exponen- 
tially with energy or some other chaoticity parameter (cf. 
the Appendix H in [s^ or, to a large extent, [13]); the 
number of NPOs is evidently even much larger. 

For the semiclassical density of states ([36|) . the sum- 
mation over POs is known not to converge in general (cf. 
jsst). For the semiclassical expressions PO)) - how- 
ever, the convergence of the sums over NPOs is appre- 
ciably improved due to the appearance of their periods 
T(A, r) in the denominators. In practice, we find that it 
is sufficient to include only a finite number of shortest 
orbits, as illustrated for example in Fig. O below. 

The expressions (gHl) - (021) are only valid if the NPOs 
going through a given point r are isolated. In systems 
with continuous symmetries, caustic points exist in which 
the Van Vleck determinant 'D± becomes singular. The 
same happens at points where bifurcations of NPOs oc- 
cur. In such cases, uniform approximations can be devel- 
oped which lead to finite semiclassical expressions: these 
will be presented in Sec. IIII E 31 and in |2]| . 

We should also emphasize that the semiclassical ap- 
proximations are not valid in regions close to the classical 
turning points r\ defined by V{r\) ~ A. Since the clas- 
sical momentum p{X,r\) in (j3ip becomes zero there, the 
spatial density (j40p always diverges at the turning points. 



Furthermore the running time T(A, r), which appears in 
the denominator of all densities PO)) - (|42p . may turn to 
zero at the turning point for certain orbits. To remedy 
these divergences, one has to resort to the technique of 
linearizing a smooth potential ^(r) around the classical 
turning points, which is familiar from WKB theory [39| . 
We shall discuss this in detail in Sec. IIVI 

Our semiclassical formulae (PU]) - (|^^ can also be ap- 
plied to billiard systems in which a particle moves freely 
inside a given domain and is ideally reflected at its bound- 
ary. The only modification is that for a given orbit, each 
reflection at the boundary contributes one extra unit to 
the Morse index fj, in (|43)) . since the difference in the 
semiclassical reflection phases between a soft and a hard 
wall is 7r/2. A detailed application of our formalism to 
the two-dimensional circular billiard, including a com- 
plete determination of all closed orbits of this system, 
has been given in [2ll |. 



C. Local virial theorem 

1. Statement and test of the theorem 

We now shall discuss a result which can be directly in- 
ferred from the semiclassical equations (|40p - (j42]) , with- 
out detailed knowledge of the NPOs that contribute to 
them in a particular potential. 

Since the modulus of the momentum p(A,r) depends 
only on position and Fermi energy, but not on the orbits, 
we have taken it outside the sum over the NPOs. Com- 
paring the prcfactors in and (|41|) and using (|31|) . we 
immediately find [l9j the relation 



5t{y) ~ [A - y(r)] 5p{Y) 



(49) 



This is exactly the local virial theorem (LVT) that was 
derived in [10| from the quantum-mechanical densities 
in the asymptotic limit iV — oo for isotropic harmonic 
oscillators. Here we obtain it explicitly from our semi- 
classical approximation. Since no further assumption 
about the potential or the contributing NPOs has been 
made, the LVT holds for arbitrary integrable or non- 
integrable systems in arbitrary dimensions with local po- 
tentials V{t) and hence also for interacting fermions in 
the mean-field approximation given by the DFT. We re- 
call, however, that p9|) is not expected to be valid close 
to the classical turning points. 

No such theorem holds for the density STi{r), since it 
depends on the relative directions of the momenta p and 
p' of each contributing orbit through the factor Q{X, r) 
P^ appearing under the sum in p2[l . 

In Fig. [2] we test P^ explicitly for the coupled two- 
dimensional quartic oscillator 

V{x,y)^^{x^ + y')-KxV, (50) 

whose classical dynamics is almost chaotic in the limits 
K = 1 and K ^ —oo pol.l4H. but in practice also for k = 
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FIG. 2: (Color online) Oscillating part of spatial densities of 
TV = 632 particles in the nearly chaotic potential (|50|l with 
K — 0.6 {h = m — 1). Top: The solid (black) line gives the 
r.h.s. of the LVT (|49|l . the dashed (red) line gives 5T{x,y), 
and the dotted (blue) line gives 5T\{x,y), all taken along the 
line y = x/^/i. Bottom: S^{x,y) along y = xjypi. 



0.6 (see, e.g., [42|)- We have computed its wavefunctions 
using the code developed in [i^. In the upper panel of 
Fig. [2] we show the left side (dashed line) and the right 
side (solid hnc) of the LVT for this system with N = 
632 particles, using the exact densities along line y = 
xj-sfZ, i.e., 6p{x,x/^f^) and 5t{x,x/^/?>). The agreement 
between both sides is seen to be very good, except in the 
surface region. We also show 5ti{x , x / ^/?>) (dotted hne). 
This demonstrates that the leading contributing NPOs 
in this system are not self-retracing. Correspondingly, 
the quantity (5^(.t, a:/\/3) in the lower panel is seen not 
to be negligible. 



D. D = 1 dimensional systems 

In a one-dimensional potential V{x) there is only linear 
motion along the x axis. As discussed in [l9| . the only 
types of NPOs are those running from a given point x to 
one of the turning points and back, including fc > full 
periodic oscillations between both turning points. We 
name the two types of orbits the orbits that start 
from any point x 7^ towards the closest turning point 
and return to x, and the "— " orbits that are first reflected 
from the farthest turning point. Clearly, these orbits have 
opposite initial and final momenta: p = —p', so that the 
momentum mismatch function (jU) equals (5(A, x) = —1. 
Consequently, one obtains from (|42)) directly the relations 



(51) 



6ti{x) :^ -5t{x) , 5£,{x)~Q. 



Note that these results do not hold near the classical 
turning points, where the semiclassical approximation 
breaks down (cf. Sec. IIVI sec also the example in Fig. 



[TOl where 5^{x) is small inside the systems but becomes 
comparable to 5p(x) near the turning points.) 

The explicit evaluation of (|40)) for D ^ \ was done 
in for smooth potentials; the result in the present 
notation is 



5p{x) 



^_^^^co^[kS'i\\x)/h] 



t'^\\x) 



(52) 



where 5j_ are the actions of the and " type NPOs 

are there running 

for the quartic 



(including k full periods), and T^'^ 
times defined by ((M]) . 

A numerical example was given in jT 
oscillator in one dimension 



V{x) = a;V4 . 



(53) 



Unfortunately, an error occurred in the drawing of Fig. 
1 in [1^; the present Fig. [3] is its corrected version. In 
the upper panel, it is seen that the semiclassical approxi- 
mation (|52[) for Sp(x) agrees very well with the quantum 
result, and in the lower panel the relations ([i^ and f?!]) 
between the quantum results are seen to be well fulfilled. 
The only sizable deviations occur very near the classical 
turning point, as expected. 




FIG. 3: (Color online) Upper panel: Oscillating part Sp{x) of 
the particle density of A''=40 particles in the quartic potential 
(|53p (without spin degeneracy; units: fi — m — 1). Dots 
(black) show the quantum-mechanical result; the solid line 
(red) shows the semi-classical result (|52p . and the dashed line 
(blue) the approximation (|65|l (for D=l) valid for small x 
values. Lower panel: Tests of relations (|49p and (|5ip between 
the quantum-mechanical densities for the same system. Solid 
line (red): St{x), dashed line (blue): —5ti{x), dotted line 
(black): r.h.s. of dMJ. [Corrected figure from [l^.] 
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We emphasize that the Friedcl oscillations near the 
surface are dominated by the primitive "+" orbit (with 
fc = 0). Its contribution diverges, however, since its run- 
ning time t!^\x,x) tends to zero there. This divergence 
can be remedied in the WKB-type linear approximation 
to the potential which we discuss in [2^ for smooth po- 
tentials, or by the short-time propagator for hard- wall 
potentials (i.e., bilhard systems) discussed in Sec. lIVBl 
First we will, however, examine the strictly linear poten- 
tial for which the WKB approximation is exact. 



2. The 1- dimensional box 

For the one-dimensional box of length L, Eq. ((52|) 
has to be modified by omitting the phase factor (— 1)''', 
since each turning point gives two units to the Morse in- 
dex. Using ptf = 2/7ry^277iXTF/fi? and summing over 
all fc, one finds that it reproduces exactly the quantum- 
mechanical p{x) in the large- iV limit, so that the semi- 
classical approximation here is asymptotically exact. 



1. The linear potential 



E. D > 1 dimensional potentials with spherical 
symmetry 



In [28[, we give the exact quantum- mechanical den- 
sities for the one-dimensional potential V{x) — ax. Al- 
though this potential does not bind any particles, its den- 
sity close to the turning point will be of use in Sec. IIV Al 
Here we give its semiclassical analysis. 

Since a particle cannot be bound in this potential, the 
only closed classical orbit starting from a point x is the 
primitive orbit "-f" (fc = 0) going to the turning point 
xx ~ X/a and back to x. Its action is 



S+ (x) = 5^ ' (x) 2 / p( A, x) dx 



3a 



= h-\zx\'/' ^h2a, (54) 

where the last equalities make use of the quantities de- 
fined as 



2m 



and 



zx = a{ax — A) , Po = 2 



1/3 



2ma 



1/3 



(55) 



2aa. (56) 



Using (j40|) for D = 1, we obtain the semiclassical contri- 
bution of this orbit to the spatial density [cf. Eq. (23) of 
[3 with a = +, k^l] 



5p{x) 



1 



27r (A - ax) 



S+{x) 



(57) 



which is identical to the asymptotic expression for the ex- 
act quantum-mechanical result [l^. Thus, the orbit "+" 
creates the Friedel oscillations. Using the LVT (|49| and 
Q = —1 in (j42p . we obtain immediately the expression 
for the kinetic-energy densities 



5t[x) = —5ti{x) = — — cos 
27r 



\s^{x) 



(58) 



which is identical to the asymptotic quantum result [28| . 
The expression ([57)1 diverges at the classical turning point 
Xx ■ To avoid this divergence one has to use the exact ex- 
pressions [i^ , which can be considered as the regularized 
contributions of the primitive "+" orbit near the turning 
points. 



In this section we discuss potentials in D > 1 with 
spherical symmetry, so that V{v:) = V{r) depends only 
on the radial variable r = |r|. The particle number N is 
chosen such that energy levels with angular-momentum 
degeneracy are filled so that all spatial densities, too, 
depend only on r. In such systems, the two kinds of 
oscillations discussed in Sec. Ill B 21 can always be sepa- 
rated clearly in the central region r ~ 0. Indeed, this 
behavior is explained by the fact that the angular mo- 
mentum of the orbits is conserved. Therefore, the shape 
of a closed orbit whose starting point r approaches the 
center of the potential tends to become flattened and con- 
centrated near a radial periodic orbit. Thus, close to the 
center there are only two types of non-periodic orbits: 
Firstly, the radial orbits of the same types "+" and "— " 
as discussed for the one-dimensional case, with opposite 
momenta p ~ — p', leading to the same kind of oscilla- 
tions that we know for D = \. Secondly, non-radial orbits 
which near r = have almost equal momenta p ~ p', so 
that they become nearly periodic. 

Semiclassically, the two types of radial and non-radial 
NPOs are responsible precisely for the two kinds of os- 
cillations which we described in Sec. IIIB 21 The regu- 
lar short-ranged oscillations, denoted (5rp(?') etc., can be 
attributed to the radial and "— " orbits. The long- 
ranged irregular oscillations, denoted 5n-T-p[r) etc., must 
be attributed to the non-radial NPOs: these lead to slow 
oscillations because their actions are almost independent 
of the starting point near r = 0. 

The contributions of the radial NPOs in radially sym- 
metric systems has already been anticipated in [l9| ; they 
will be discussed in the following section. In particular, 
like for D = 1, the primitive orbit is seen to be solely 
responsible for the Friedel oscillations near the surface of 
a, D > 1 dimensional spherical system. 

Non- Radial orbits can only occur if there exist classical 
trajectories which intersect themselves in a given point 
r. As is well known from classical mechanics, such or- 
bits do not exist in isotropic harmonic oscillators (and in 
the Coulomb potential). This explains the fact that no 
irregular long-ranged oscillations are found in the den- 
sities of harmonic oscillators (or, trivially, in any 
onc-dimensional potential) . 

We emphasize that for D — 2, all closed NPOs are 
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isolated except if they start at r = 0, in which case they 
form degenerate famihes due to the radial symmetry (cf. 
Sec. IIII E 3| . In I? > 2 dimensions, however, also the 
non-radial NPOs starting at r > have continuous ro- 
tational degeneracies. For the corresponding families of 
orbits, the Van Vleck determinant T)-^ in the semiclas- 
sical Green function (P^)) becomes sin gula r at all points 
r. This divergence can be removed |43| by going one 
step back in the derivation of (p9)) . In the convolution 
integral for the time-dependent propagator, one has to 
perform a sufficient number of intermediate integrals ex- 
actly rather than in the stationary-phase approximation 
(for details, see [i^ where this was done to obtain the 
trace formula ([36]) for systems with continuous symme- 
tries). As a result, the semiclassical amplitudes of the 
degenerate families of orbits are of lower order in h than 
for isolated orbits and thus have a larger weight. In our 
present case, the h dependence of the ratio of amplitudes 
between the irregular and the regular oscillations e.g. in 
the particle density becomes: 



l'^irrP(?-)| 

\5.p{r)\ 



cx fl 



{D>1) 



(59) 



The same ratio holds also for the other spatial densities. 
This can be seen, e.g., in Fig.[T]for the spherical billiard in 
-D = 3, where the amplitude of the irregular oscillations 
is larger than that of the radial oscillations (except near 
r = 0). In passing, we note that for spherical billiards 
with radius R, the energy dependence of the semiclassical 
results scales with the dimensionless variable p\R/h, and 
the ratio ^ becomes |i5irrp/^rp| oc [pxR/h]'^^-^'^/^ . 

It should be stressed that the separation of two classes 
of NPOs and hence the two types of oscillations is not 
possible in systems in D > 1 dimensions without radial 
symmetry. This will be illustrated in Sec. IIIIFI 

A further complication in systems with Z? > 1 is that 
the NPOs can undergo bifurcations under the variation 
of the starting point r. At these bifurcations, new NPOs 
or POs are created. This is discussed extensively in a 
publication [2lj on the two-dimensional circular billiard. 
For this system, a complete classification of all NPOs 
could be made and analytical expressions for their actions 
and Van Vleck determinants have been derived. 



1. Contributions of the radial orbits: Earlier results 

Recall that since all radial NPOs fulfill p' = — p, they 
have Q{X,r) ~ —1 under the sum in (|42p. Therefore we 
immediately obtain the semiclassical relation 



(60) 



Indeed, this was found to be fulfilled, sufficiently far from 
the turning point, for all quantum systems discussed in 
SccHlBl 

In order to derive some of the other forms of local virial 
theorems discussed in Sec. IIIBl it is important to notice 



the action of the differential operator V on the semiclas- 
sical density in ([^0]) . The contributions of leading-order 
in h (i.e., the terms of the largest negative power of h) 
come from the phase <I>(A, r) given in (|43l) . From the 
canonical relations we find 



yg«*(A,r) ^ i (p' _ p) e'*(^.'-) 



(61) 



y2 g.*(A,r) ^ _ J_ _ p)2 g^$(A,r)^ (g2) 



and 



which occurs for each NPO under the summation in (|40p . 
For the radial orbits, one therefore obtains with (|3ip the 
following differential equation for 5^p{r), which was al- 
ready given in p^ : 

-^VH,p{r)c^[\-V{r)]5,p{r). (63) 

For small distances r from the center so that V{r) ^ A, 
(I63p becomes the universal Laplace equation 



(64) 



which was obtained asymptotically from the quantum- 
mechanical densities of isotropic harmonic oscillators in 
p^ . It has the general solution 



5rP{r) = (-1) 



Ms- 1 m 



hT,i{X) 



{£^JjA2rpx/n). (65) 



Here Jij{z) is a Bessel function with index i/ — D/2 — 1, 
Ms = + 1 is the number of filled main shells [i^l, Tli is 
the period of the primitive radial full oscillation and px = 
(2mA) "'^/^ is the Fermi momentum. The normalization of 
cannot be obtained from the linear equation (|64p : 
we have determined it from the calculation presented in 
Sec. IIII iTSl For harmonic oscillators, where T,! = 2tt/uj, 
equation (p5|) becomes identical with the result in \W^ . 
Eq. (69), that was derived from quantum mechanics in 
the large- iV limit. 

The quantity 5,:p{r) can also be calculated directly 
from (|40p . including only the radial NPOs. The summa- 
tion over their repetitions goes exactly like in the one- 
dimensional case done in [l9|, except for the evaluation 
of the determinant 7)^. This determinant becomes sin- 
gular at r = due to the continuous degeneracy of the 
and "— " orbits: the point r = is a caustic point for 
all radially symmetric systems with D > 1. The regular- 
ization of this singularity, leading precisely to the result 
jS]) . is discussed in Sec. IIII E3] below. 



2. Isotropic harmonic oscillators in D dimensions 

We now investigate the densities in the isotropic har- 
monic oscillator (IHO) potential in D dimensions defined 
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V{r) 



,2^2 



r e 



(66) 



First wc mention the well-know fact that in IHO poten- 
tials with arbitrary Z? > 1, all orbits with nonzero angu- 
lar momentum are periodic, forming ellipses which may 
degenerate to circles or radial librations. Hence the only 
NPOs are the radial orbits and "— " . Since we just 
have seen that in the leading-order semiclassical approx- 
imation, (5rTi(r) = —5i-T{r), it follows that (5^(r) = to 
leading order like for D = 1, thus explaining the smooth 
behavior of for IHOs 

For the IHO potentials, the transverse determinant Pj^ 
can be easily computed. It is diagonal and reads 



\V±{X,r)\ 



mX 



rp{\, r) 



D~l 



(67) 



which does not depend on the type and the repetition 
number k of the orbit. Following (|40)l and [l^, we com- 
pute Sp{r) as a sum over the contributions of the 
and "— " orbits, which is given by 



dp{r) 



1 



(27rft)^ p(A,r) 



mX 



rp{X, r) 



<E 

k=0,± 



COS 



si'\x,r)~{D + l)^~f^^^'>^ 



Ti'\x,r) 



(68) 



Here we have used the analytical form of the actions and 
periods 




FIG. 4: (Color online) Oscillating part of the spatial particle 
density times for AD IHO for = 632502, i.e. with M = 50 
filled shells (units: fi = m — uj — 1). Dots are the quantum 
results. The solid (red) line is the analytical expression (|68p 
using the Morse indices given in (|7ip . and the dashed (blue) 
line is the asymptotic formula (|65p valid close to r = 0. 



is a caustic point due to the spherical symmetry. This 
divergence will be regularized in the following section. 

Using the Morse indices ([7T|) and for A the expression 
X = hu! [M + {D + l)/2] [l3|, we can perform the sum- 
mation over k in (j68p analytically for small r, like it was 
done in [l^ for the ID case. The result then is exactly 
that given in ([65|) with Tri(A) = 27r/ijj, but replacing the 
Bessel function Ju{z) by its asymptotic expression for 
large argument z, i.e., using 



Ju{z) 



C0S(Z — VTT /2 — TT / A) . 



(72) 



S'i\x,r) 



T'i\x,r) 



ttX ~ 
(2fc+ 1) — Trp{X,r) 

UJ 



2X . ( mu)r 
=p — arcsm 



UJ 



, N 2 . / rriLor 
{2k + 1)— =F — arcsm 

UJ LU \ P\ 



(69) 
(70) 



We compute the Morse indices following Gutzwiller [4a |. 
Each turning point contributes a phase of 7r/2. Besides 
we evaluate the number of extra conjugate points includ- 
ing their multiplicities depending on the dimension, con- 
tributing a phase tt{D — l)/2 each (they are most easily 
determined from the propagator of the harmonic oscilla- 
tor in the time representation). The final result for the 
Morse indices is 



(k) 



2kD + 1 , 



(k) 



2kD + D . 



(71) 



We note that the equation ((68|) is consistent with re- 
sults derived in [ll[ from the quantum mechanical den- 
sity p{r). Fig. [4] shows a comparison of the semiclassical 
results (|55)) with the exact quantum result for the case 
D — 4:. We have multiplied both by a factor since the 
semiclassical determinant diverges at r = which 



3. Regularization close to the center 

In this section we compute the contribution of radial 
NPOs to the semiclassical particle density close to the 
center of an arbitrary potential with radial symmetry. 
As stressed in the last section, the semiclassical Green 
function for D > 1 is not defined at r = where 'D± di- 
verges. The reason is the caustic that occurs there: fixing 
the position of the point r = r' = does not uniquely 
determine a closed orbit (periodic or non-periodic) which 
belongs to a continuously degenerate family due to the 
spherical symmetry. A standard method to solve this 
problem is to introduce the mixed phase-space represen- 
tation of the Green function close to the diverging point, 
as proposed initially by Maslov and Fedoriuk [46| . 

Here we follow more specifically the procedure outlined 
in [47j . The mixed representation of the Green function 
can be approximated in a form analogous to that in the 
coordinate representation. This is due to the smoothness 
of the phase-space torus which implies that no diverging 
points can occur simultaneously in position and momen- 
tum (cf. [i^). Following Gutzwiller, we use for every 
classical trajectory 7 an "intrinsic" (or local) coordinate 
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system r = (r||,rj^). where the coordinate r|| is taken 
along the trajectory and is the vector of all other 
coordinates transverse to it; p = {p\\,p±) is the corre- 
sponding system for the momentum. We next re-write 
the coordinate-representation of the Green function as 
the inverse Fourier transform of the mixed Green func- 
tion with respect to the final transverse momentum p'^: 

G,ei(£;,r,rf|,rl) = ^ J dp'^ 

X G^{E, r, r[| , pi) cxp (^^ r'^ • pi) , (73) 

where the sum is over all classical trajectories 7 starting 
at r and ending at {rpp'_^_) in phase space. Hereby the 
contribution of the orbit 7 to the semiclassical mixed 
representation of the Green function is given by [i^: 



G^{E,r,r\^,p'j_) = aDV^{E,r,r\^,p'j^) 
X exp ( ^Sj{E, r, rj, , pj 



^^^{E,r,^,p^) - —n^ 



(74) 



where S is the Legendre transform of the action S be- 
tween the variables and p^: 



S^{E, r, rj, , p1) = S^{E, r, r\ , r^) - v\ ■ p' 



(75) 



Since in the mixed-representation Green function, we 
have to evaluate the action S for radial orbits with fixed 
momentum close to the center, the rotational symmetry 
in position is removed and G is regular. The Van Vlcck 
determinant in this representation is 



1/2 



biip'il^/' 

and the Morse index becomes 



2?j_-y = det 



dpi_ 
dp' 



(76) 



firy for positive eigenvalue of det 
yU-y + 1 for negative eigenvalue of det 



\dp'^ 



dp\ 



Far from singular points, the evaluation of (|73|| using 
the stationary phase approximation (SPA) yields [4g the 
standard semiclassical Green function (P^ . 

After performing the h expansion and the integration 
over the energy similarly as in (Toj . the oscillating part 
of the particle density is given by: 

5p{v) = 2^Im|^ j Ap\G{X,v,^,p\) 



(77) 



X exp [ -rj_ • p 



Close to the center of the potential, now, we replace the 
non-radial NPOs by the radial ones with 7 ="±" orbits 



with fc-th repetitions. For non-periodic orbits in the ra- 
dial direction r we have r|| = r. We neglect the higher 
orders in ri, leading to the following approximations: 

|det(9pi/9pV)| « 1, 

« P||(A,pV):= ^2m\-p'l , 
(r,ri) « (r,0). 



S. 



ik) 



(k + 1/2) 5ri T 2rp' 



T. 



« (fc + l/2)r,: 



(78) 



Furthermore, we approximate the action S^i of the prim- 
itive periodic diameter orbit by S^i w 27rfi,[A/ + (D + 
l)/2]. This is exact for IHOs where S-ci = 27rA/a; and 
A = fttj [M + (D + l)/2] can be used [l3|; for arbitrary 
radial potentials it corresponds to a radial WKB quan- 
tization, whereby M is a "main shell" quantum number 
that has to be suitably chosen [31 ■ Also, we assume that 
each eigenvalue of det((?r'x/9p'^) is negative (positive) 

for the orbits ("—"), leading to /i^'' = /i^''. This is 
again exact for IHOs; for other radial potentials we have 
verified its validity numerically. With these approxima- 
tions, the sum over the repetitions of all radial orbits can 
be performed exactly like in the previous section. The 
oscillating part of the particle density then simplifies to: 



(-1) 



dp 



, co s[2rp|[(A,p'j_)/^] 

P||(^:P'±) 



(79) 

The integration has to be taken over half the solid angle 
in the {D — l)-dimensional transverse momentum space, 
avoiding a double-counting of the two orbits. So it is 
natural to make a change of variables to dimensionless 
hyper-spherical coordinates. Using the integral represen- 
tation of the Bessel functions 14911 



Ju{z) 



2(z/2)- 



cos(zt)dt, (80) 



we obtain exactly the same result as in (j65p . confirming 
its normalization. 

We stress that this regularization is only valid near 
the center, i.e., for r ~ 0, as can be seen in the exam- 
ple of Fig. m where the result is displayed by the 
dashed line. The reason is that for larger values of r, the 
approximations (j78p are no longer valid. If one restricts 
oneself to the leading contributions of the primitive orbits 
and "— " with fc = 0, a "global uniform" approxima- 
tion can be made which interpolates smoothly between 
the regularized result ([S5)) near r = and the correct 
semiclassical contributions obtained from (|40l) at larger 
r. This uniform approximation is derived and used in 
pll | for the 2D circular billiard system which we briefly 
discuss in the following section. 
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4- The two-dimensional circular billiard 



The two-dimensional circular billiard, which can be 
taken as a realistic model for quantum dots with a large 
number N of particles, has been investigated semiclassi- 
cally in (2l| , where all its periodic and nonperiodic closed 
orbits have been classified analytically. We discuss there 
also the various bifurcations at specific values of the ra- 
dial variable r, at which POs bifurcate from NPOs or 
pairs of NPOs are born. At these bifurcations, the semi- 
classical amplitudes in - must be regularized 
by suitable uniform approximations. We refer to [2l| for 
the details and reproduce here some numerical results to 
illustrate the quality of the semiclassical approximation. 
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FIG. 5: Particle density in the two-dimensional disk billiard 
with radius R, containing — 606 particles (units: h^/2m — 
R= 1), divided by TV. The solid line is the quantum result, 
the dotted line the semiclassical result with all regularizations 
(see [U for details). 



Fig.[5]shows the total particle density p{r) for N ~ 606 
particles in the circular billiard. The solid line gives the 
quantum result, obtained from ^ using the solutions of 
the Schrodinger equation with Dirichlet boundary con- 
ditions, which are given in terms of cylindrical Bessel 
functions. The dotted line gives the semiclassical re- 
sult, obtained by summing over the ^ 30 shortest NPOs. 
(Hereby we used the regularization of the radial "-f " and 
"-" orbits at r = by that of the primitive "+" 

orbit near r = i? by ([55)1 given in Sec. IIVB II below, 
and uniform approximations for the bifurcations of some 
of the non-radial NPOs as described in detail in [2l[.) 
We see that, indeed, a satisfactory approximation of the 
quantum density can be obtained in terms of the shortest 
classical orbits of this system. 

In Fig. [6] we demonstrate explicitly the contributions 
of non-radial NPOs to the kinetic-energy densities ri(r) 
and ^(r) close to the center, calculated as in Fig. [5] but 
for N — 9834 particles. We clearly see that <5^(r) is not 



FIG. 6: Oscillating parts of kinetic-energy densities, (5ri(r) 
(fast oscillations) and 5£,{r) (slow oscillations) for A'' = 9834, 
divided by N^^^ . Solid lines: exact quantum results. Dashed 
lines: semiclassical results and units as in Fig. (5) 



smooth; its slow, irregular oscillations arc due to non- 
radial NPOs which have the form of polygons with 2k 
reflections (fc = 1, 2, . . . ) at the boundary and one corner 
at a point r close to the center. The first fcmax = 20 of 
them were included with the appropriate regularization 
at r = where they are degenerate with the fc-th repeti- 
tions of the diagonal PO (see [2l| for details). The agree- 
ment between quantum and semiclassical results is again 
satisfactory; the discrepancy that sets on for r > 0.18 is 
due to the missing of more complicated non-radial orbits. 
The quantity (5Ti(r), on the other hand, clearly exhibits 
both kinds of oscillations according to (PS)) : the slow ir- 
regular part, which is identical with (5^(r), is modulated 
by the regular fast oscillations due to the radial orbits. 



F. D > 1 dimensional systems without continuous 
symmetries 

In D > \ dimensional systems without continuous 
symmetries, it is in general not possible to find the clas- 
sical orbits analytically. As in POT, the search of closed 
orbits must then be done numerically. A practical prob- 
lem in such systems is also that the densities as func- 
tions of D coordinates are not easily displayed. For tests 
and comparisons of various approximations or of the lo- 
cal virial theorems, we have to resort to taking suitable 
one-dimensional cuts (i.e., projections) of the densities. 
In the following paragraph, we discuss a class of inte- 
grable billiard systems, in which all closed classical orbits 
can easily be found and their semiclassical contributions 
to the densities can be analytically obtained. These are 
ZJ-dimensional polygonal billiards that tessellate the full 
space under repeated reflections at all borders. We illus- 
trate the method for the example of a rectangular bil- 
liard. Although this does not correspond to any physical 
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system (unless experimentally manufactured as a rectan- 
gular quantum dot with many electrons), it is a useful 
model without spherical symmetry that allows for ana- 
lytical calculation of the classical orbits and their prop- 
erties. 



For billiards, classical trajectories are straight lines 
which are reflected at the boundary according to the 
specular law. Let us consider a two-dimensional billiard 
that tessellates the plane, such as the rectangular bil- 
liard shown in Fig. [71 Choose a trajectory starting at 
a point P, reflected at the point Rq and reaching the 
point Pi. Now reflect the boundary at a side containing 
the point Rq, the image RoP{ of the segment RqPi gives 
the straight line PP{. The next portion of the trajec- 
tory after reflexion in Ri , can be found by reflecting the 
new billiard at the side containing R'l. This process can 
be repeated until the trajectory ends. To get the closed 
trajectories at P we have to compute all images of P in 
the images of the billiard. Now a straight line joining P 
and an image of P gives a closed orbit. Thus, construct- 
ing all images of P by simple geometry, enables one to 
compute all trajectories and their related initial and final 
momenta for I?-dimensional polygonal billiard that fills 
the D-dimensional Euclidian space. Note that the Jaco- 
bian for these systems is easily computed and equals 
(p/iNPo)^~^ where Lnpo is the length of the orbit. 



FIG. 7; Images (triangles and crosses) of a point P{x,y) (full 
circle) for a rectangular billiard. Joining by a straight line the 
full circle to a cross gives a non-periodic orbit whereas joining 
to a triangle gives a periodic orbit. 

We illustrate this method for the case of a 2D rectan- 



gular billiard with side lengths Qx and Qy. There are 
four types of images of P{x,y); one leading to POs and 
three (labeled by the index a, b and c) leading to NPOs 
(see Fig [7]). Table U lists the basic ingredients to compute 
the spatial densities, using L{x,y) = 2y^ + and 



1. Billiards tessellating flat space: the rectangular billiard fi^TUjI^) — 



■ cos 



PxL{x,y) Stt 



(81) 

with PA = (2mA)i/2. From (gO]), (gS]) for D = 2 we obtain 

OO 

5p{x,y) = ^Pi(^'y)^ (82) 

kj:,ky— — oo /— a,b,c 
oc 

STi{x,y) = X! X! ^■^i'(^'y)' (83) 



^— a.b.( 



where the partial contributions 5pi{x,y) and 6Tii{x,y) 
for the orbits of types / = a, b, and c are given in Tab. 
m 5t{x, y) is obtained from ([5^ using the LVT 
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FIG. 8: (Color online) Oscillating part of the spatial densities 
for a rectangular billiard with sides Qx ~ 2^'''' and Qy = 3^'''* 
for N = 2000 along the line x = Qx/2 (units: ft^ = 2m = 
1). Solid (black) lines are the semiclassical results using (|82l) . 
(|83p and ((49]); dashed (red) lines are the quantum- mechanical 
results. 

We now present numerical results for the rectangular 
billiard with side lengths = 2-*^/^, Qy = 3^^^ (units: 
— 2m ~ 1), containing = 2000 particles. In Fig. 
Owe show the quantities Sp (top), St (center) and 5ti 
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y 

FIG. 9: (Color online) Same system as in Fig. [S] Here, se- 
lected contributions to (|82p of the primitive NPOs with fc = 
are shown. The full (black) line gives the contributions of the 
primitive self-retracing orbits, and the dashed (red) line that 
of all other primitive orbits. 



(bottom) as functions of y with fixed x = Qx/'^- Dashed 
lines are the quantum- mechanical results, solid lines the 
semiclassical ones using (f82|) . (f83l) and (|49|) . We see that 
summing over all orbits yields very good agreement, ex- 
cept close to the boundary where the Friedel oscillations 
were not regularized. 

In Fig. [9] wc display selected contributions of some 
of the primitive orbits (fc = 0) to the particle density 
Sp{x,y). The solid line gives the contribution of self- 
retracing orbits with p = — p', and the dashed line that of 
the other primitive NPOs. It is evident that no clear sep- 
aration of regular short-ranged and irregular long-ranged 
oscillations can be made here. 



IV. REGULARIZATION NEAR SURFACE 

As we have pointed out in the previous section, the 
semiclassical approximation of density oscillations in 
terms of classical orbits breaks down near the classical 
turning point due to the diverging amplitude of the prim- 
itive orbit (with k = 0) which close to the surface is 
responsible for the Friedel oscillations. In order to regu- 
larize this diverging amplitude, different techniques must 
be used for smooth potentials and for billiards with re- 
flecting walls. 



A. Smooth potentials 

In smooth potentials ^(r), the divergence can be reg- 
ularized by linearizing the potential near the classical 
turning points, as it is done in the standard WKB ap- 
proximation (39[ . In the surface region close to a tur ning 
point, the exact results for linear potentials given in |28| 
can then be used. We demonstrate this first for the one- 
dimensional case, and then illustrate it also for potentials 
in 13 = 3 with spherical symmetry. 



1. Linear approximation to a smooth ID potential 

We start from an arbitrary smooth binding potential 
V{x) and approximate it linearly around the turning 
point x\ defined by V{x\) ~ X. Without loss of gen- 
erality, we assume x\ > 0. Expanding V{x) around xx 
up to first order in x — x\, we get the approximated po- 
tential 

V{x) = X + a{x-xx), a = V'{xx)>0. (84) 

We can therefore apply the results of [11] . The oscillating 
part of the density near the turning point then becomes 



Spiin{x) = Po\ [Ai'(zA)]^ - 2:AAi^(zA) 



— V-zx (^{xx - x) 

TT 



where the last term is the subtracted TF part and 



Po 



zx = —\x- Xx) ■ 



(85) 



(86) 



The oscillating parts of the kinctic-encrgy densities t(x) 
and ^ix) become in the same approximation 



5riin(a;) 



2a 



Ai{zx)Ai'{zx) - zx[Ai'{zx)]' 



+zlAi'{zx)--\zx\'^'eixx-x) 



(87) 



Aiizx)Ai'{zx) + 2zx[Ai'izx)]^ 
-2z2Ai2(zA) + -|2Ap/'e(xA-x) 



In the next step, we introduce uniform linearized approx- 
imations, in which the argument zx in ((85l) . ((87)) . and ((88] 
is not as given in but replaced by 



Zx 



[3S+ix)/4hf/^ , 



(89) 



where S+{x) is the correct action of the "-|-" orbit for 
the given potential V{x). On can show that this relation 
is exact for the linear potential; it is uniform for other 
smooth potentials in that it holds locally at the turning 
point and yields the correct phase of the oscillation at all 
other distances from the turning point. 

Figure[Tn]shows numerical results for these uniform ap- 
proximations for the quartic oscillator ([53]) with iV = 40 
particles, compared to the exact quantum results. We 
see that the uniform linearized approximation reproduces 
very well the Friedel oscillations near the turning point in 
all three densities. The phase of the oscillations is seen to 
be correct at all distances. The amplitudes are not exact 
in the asymptotic region, i.e., near x = 0. This is not 
surprising, since the contributions of all "— " orbits and 
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FIG. 10: Oscillating parts of densities for the quartic potential 
(|53|) with N = 40 particles (units h = m — 1), shown on the 
same scale. Solid lines: exact quantum-mechanical results, 
dotted lines: uniform linearized approximations (|85p . (|87p and 
(f88)l with the argument zx given in (f89)l. 

those of the "+" with fc > are missing in this approx- 
imation. We see that S^{x) vanishes inside the system, 
as expected from the semiclassical leading-order resuh 
on the r.h.s. of (|5ip . However, near the turning point, 
where the semiclassical approximation breaks down, the 
magnitude ofd^{x) is comparable to - and for the quartic 
potential even larger than - that of 6p{x). (Note that all 
three density oscillations are shown on the same vertical 
scale.) 

2. Linear approximation to smooth radially symmetric 
potentials in D > 1 

We now start from an arbitrary smooth binding po- 
tential with radial symmetry, V{r) = V{r), r = |r|, in 
D > I dimensions. As above, we replace it by its linear 
approximation around the turning point rx analogously 
to dMl): 

F(r) = A + a- (r-TA), a = VV{rx). (90) 

Due to the spherical symmetry of l^(r), all components 
of the vector a have the same magnitudes: 

a = arx/rx, a^V'{rx). (91) 

Therefore, we may choose the radial variable r along any 
of the Cartesian axes Xi , and the results 

P(^^) = -l^Pio\^Kzi\)^i'{zix) + 2zix[Ai'{zix)]^ 
4<37r 

-24Ai2(^,A)|, {D = 3) (92) 

where piQ = 2cFiai, taken from psj for the linear potential 
with D > \ apply with the replacements Xi r, Zix — > 



a{ar — A). As in Sec. IIV A II for D = 1, we may then 
subtract their ETF contribution. Finally we introduce 
the uniform approximation to their oscillating parts near 
the surface with the argument (j89p expressed in terms of 
the action S+{r) of the primitive radial orbit of the 
given radial potential V{r): 

zx = ~[3S+ir)/Ahf\ (93) 
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FIG. 11: Oscillating part of particle density for the 3D IHO 
with N = 22960 particles {Ms = 40) (units h = m = u = 1). 
Solid lines: exact results, rfas/ieci lines: uniform linearized 
approximation (|92|) from [2a | with argument 2a given in (|93|l . 
Upper panel: smooth part in Sp{r) taken as TF density, lower 
panel: smooth part in Sp{r) taken as ETF density. 

In Fig. [TT] we show numerical results for this approx- 
imation for the 3-dimensional IHO with Ms = 40 occu- 
pied shells. The upper panel shows the exact result for 
Sp{r) (solid line), whereby only the TF approximation 
was used for its smooth part: Sp{r) = p{r) — ptf('')- 
We notice that the oscillations in the interior are not 
symmetric about the zero line, which is due to smooth 
errors in the TF density. In the lower panel, the ETF 
corrections have been included in dp(r); now the oscil- 
lations are symmetric about zero. The price paid for 
this is that Sp{r) diverges at the classical turning point. 
The uniform linear approximation (|92p with the argu- 
ment shown in both panels by the dashed lines, 
reproduces well the Friedel oscillation near the surface. 
In the interior, it fails due to the missing contributions 
of the repetitions (fc > 0) of the "-f" and of all "— " 
orbits. Once more, these results demonstrate that the 
Friedel oscillations near the surface are semiclassically 
explained by the primitive orbit alone. Its diverging 
amplitudes according to (|40p must, however, by regular- 
ized by the uniform linear approximation. In Fig. 1121 we 
show the total density for the 3D IHO with Ms = 20 
filled shells. The solid line is the exact quantum result 
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FIG. 12: Total particle density for the 3D IHO with N = 3080 
particles (Ms = 20) (units h = m = uj = 1). Solid lines: 
exact result. Crosses: semiclassical result for Sp{r) in (|68[l . 
summed up to fcmax = 15, plus petf(?'). Dashed line: uniform 
linearized approximation (|92p from [281 1 with argument in 
IMl). 



(x) is not a simple oscillatory function of the energy. 
In the following, we give results for the contributions to 
the particle density p(r) in two special cases. Unfortu- 
nately, we have not been able to derive the corresponding 
contributions to the kinetic-energy densities. 



1. Arbitrary 2D billiard 

For billiards in D~2 dimensions with arbitrary bound- 
aries, the uniform contribution to the particle density 
becomes (see [S^l for details): 



p\Ji{2dpx/h) 



2iihd^l - d/R 



(95) 



where d is the distance from the boundary and R its cur- 
vature radius at the reflection point. Hereby it is assumed 
that d is small enough so that there is only one or- 
bit going to the boundary and back to the given starting 
point. Note that the curvature radius R is negative if the 
boundary is convex at the turning point. 



(21). The crosses give the semiclassical result as the sum 
p-ETvir) + Spir), where the latter is calculated from the 
sum over the NPOs in up to fcmax = 15. We see 
that the semiclassical result reproduces very accurately 
the exact result up to r ~ 5.9, which is rather close to 
the turning point r\ ~ 6.48 where it diverges. The lin- 
earized approximation is shown by the dashed line; it 
approximates the exact density closely above r ^ 5.8. 
Thus, switching from the semiclassical approximation to 
the linearized one around r ~ 5.85 allows one to obtain 
a very good approximation of the density in all points. 



2. Spherical billiards in D dimensions 

For spherical billiards in D dimensions with radius R, 
the energy integral over (|94p can also be performed, and 
the regularized contribution of the primitive "-f " orbit 
becomes 



where p^ is the TF density given in (|15p . and 

z^2{R-r)px/h. 



u-l/2 



(96) 



(97) 



B. Billiard systems 

In billiards with reflecting walls, the above lineariza- 
tion is not possible since the slope of the potential is 
always infinite at the classical turning points. The am- 
plitude of the primitive "+" orbit can in such systems be 
regularized by using the following uniform approximation 
of the Green function for short times [31], ^01 : 



ih{2nh)D/2 



X H^^j^_,i^S/h-pn/2j, (94) 

where hI^'' (x) is the Hankel function of the first kind. 
To evaluate the corresponding uniform approximation for 
the particle density, we have to take the imaginary part 
of (|94p and perform the integration over the energy. This 
last step is not easily done analytically in general, since 



E 



-^det.^P-^ 



dr' 



1/2 



For D = 3, the expression ([96]) agrees with a result 
derived by Bonche [5l| using the multiple-reflection ex- 
pansion of the Green function introduced by Balian and 
Bloch [52|. For D = 1 (one-dimensional box), the result 
([M)) is also found from the exact solution. 

As mentioned above, the contribution (|96p is respon- 
sible for the Friedel oscillations in the densities near the 
boundary r = R. It is interesting to perform the spatial 
integral of ([Ml) over the volume of the billiard. Using the 
formula ([53|, 6.561.14, with p = —v) 



'-^dx^l- 



r(i/2) 



2^^ T{iy + 1/2) ' 



(98) 



the integral can be done in the limit p\ oo (i.e., 
for large particle numbers), and the asymptotically lead- 
ing term yields the following contribution to the particle 
number: 



6Ns 



1 



2j,D/2f^D-l Y{D) 



Pa , (99) 
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where Sd is the hypersurface of the D-dimcnsional 
sphere: 



Sn = 



2tt^ 
r{D/2) 



R 



D-l 



(100) 



We note that ([99|) corresponds precisely to the surface 
term in the Weyl expansion [s^ l of the particle number 
A^. The Fermi energy Atf in (fT5|) hereby has to be re- 
placed by the corresponding quantity Awcyi obtained by 
integrating the Weyl-expanded density of states to the 
particle number N . 

The role of the "+" orbit in contributing the surface 
term to the Weyl expansion of the density of states has 
been demonstrated by Zheng [5^ for arbitrary billiards 
in D = 2 dimensions. 



V. GENERAL RESULTS FOR FINITE 
FERMION SYSTEMS 

After having presented our semiclassical theory for spa- 
tial density oscillations and tested it in various model po- 
tentials, we shall now discuss some of its results in the 
general context of finite fermion systems. Besides the 
trapped fermionic gases [l| mentioned already in the in- 
troduction, we have in mind also self-bound molecular 
systems with local pseudopotentials, such as clusters of 
alkali metals (56j . treated in the mean-field approach of 
DFT with the local density approximation (LDA), for 
which the KS-potential is local [lg| . 



A. Local virial theorem 

One of our central results was given in Eq. (|49p which 
we repeat for the present discussion: 



ST{r) ~ [A - V{r)] Sp{r) . 



(101) 



We call it the "local virial theorem" (LVT) because it 
connects the oscillating parts of the kinetic and potential 
energy densities locally at any given point r. While the 
well-known virial theorem relates, both classically and 
quantum- mechanically, integrated (i.e., averaged) kinetic 
and potential energies to each other, the LVT in (jlOip 
does this locally at any point r. We recall that hereby 
the Fermi energy A of the averaged system is defined by 
Eq. (gSD. 

Since no particular assumptions need be made [57t to 
derive poip semiclassically from the basic equations ([^Dl) 
and (|4T|) . the LVT holds for arbitrary local potentials, 
and hence also for systems of interacting fermions in the 
mean-field approximation given by the DFT-LDA-KS ap- 
proach. This is in itself a interesting basic result. It may 
also be of practical interest, because it allows one to de- 
termine kinetic energy densities from the knowledge of 
particle densities that in general are more easy to mea- 
sure experimentally. We leave it as a challenge to the 



condensed matter community, in particular those work- 
ing with trapped ultracold fermionic atoms, to verify the 
LVT experimentally. 

Other forms of local virial theorems have been derived 
in [l(| from the exact quantum-mechanical densities of 
isotropic harmonic oscillators in arbitrary dimensions. 
A Schrodinger-like (integro-)differential equation for the 
particle density p{r) has also been derived in [l^. It 
would lead beyond the scope of the present paper to dis- 
cuss these results and their generalization to arbitrary lo- 
cal potentials based upon our semiclassical theory. This 
will be done in a forthcoming publication [2^ , where we 
also give exact expressions for spatial densities in linear 
potentials of which we already have made use in Sec. IIVI 



B. Extended validity of the TF kinetic-energy 
functional 



Presently we discuss the direct functional relation (fT^ 
between the particle and kinetic energy densities ob- 
tained in the Thomas-Fermi model. While Eq. is 
exact only when applied to the TF expressions ^T5\i and 
([T7]) of the (smooth) densities, we shall now show that a 
semiclassically approximate relation holds also between 
the oscillating exact densities: 



r(r) ~ TTF[p(r) 



(102) 



Eq. (|102p states that the TF relation IH]) holds approx- 
imately, for arbitrary local potentials ^(r), also for the 
exact quantum-mechanical densities including their quan- 
tum oscillations. This had been observed numerically 
already earlier [5^, but without understanding of the 
reason for its validity. 

The proof of (|102|) is actually very easy, having the 
LVT (|101[) at hand. Inserting p(r) = pTF(r) + Sp{r) into 
(|17p and Taylor expanding around /9TF(r), we obtain 



TTF[p(r)] = TTF[PTF(r)] 



d TTF [p] 



dp 



-0[{Spf 



Spir) - 

PTF(r) 

(103) 

Using the obvious identity ttf [ptf (r)] = ttf (r) and the 
fact that dTTF[pTF(r)]/dpTF(r) = [A — ^^(r)], we see im- 
mediately with (|10ip that, to first order in the oscillating 
parts, we have indeed the relation 



rTF[p(r)] ^ TTF(r) + ST{r) 



(104) 



We stress that, although the TF expression for all three 
kinetic energy densities T(r), Ti(r) and ^(r) is the same 
[cf. mi)], the relation (fT02)) holds only for r(r). The rea- 
son is that the LVT also only holds for this kinetic energy 
density, as discussed explicitly in the previous sections. 

In Figs. [12] and [T3] we present numerical tests of the 
relation p02[) for the two-dimensional coupled quartic 
oscillator (|50p . which represents a classically chaotic sys- 
tem, with two different particle numbers. An example for 
the three-dimensional spherical billiard, which is a good 



18 



approximation for the self-consistent mean field of very 
large alkali metal clusters [ll], is shown in Fig. [T5l We 
see that in all cases, the relation (|102p between the exact 
quantum-mechanical densities r(r) and p{r) is extremely 
well fulfilled; only close to the classical turning points, 
where the LVT (jlOip does not hold, do we see a slight 
deviation. Obviously, the terms of order O [(^p)^] , ne- 
glected in the above derivation, play practically no signif- 
icant role in the interior of the systems - even for mod- 
erate particle numbers N as seen in Fig. [T3] or in the 
examples given in Ref. (and reproduced in Q). 
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FIG. 13: (Color online) TF functional relation (fT02ll for the 
same system as in Fig. [2](Af = 632 particles). Cuts along the 
diagonal x — y. The solid (black) line is the l.h.s., and the 
dashed (red) line is the r.h.s. of (|102p . 




0.5 1 1.5 2 2.5 3 

X 



FIG. 14: (Color online) Same as in Fig. [13] for iV = 42 parti- 
cles. Cuts along the line y = x/VS. 

This result might come as a surprise, since it is well 
known from the ETF model that for smooth densities the 
gradient corrections to the functional ttf[p] do play an 
important role for obtaining the correct average kinetic 
energy. (For three-dimensional systerns, the first of them 
is the famous Weizsacker correction [531 ■) Examples for 
this are given in chapter 4.4 of [1^. However, if gradient 
corrections up to a given order were consistently added to 
()102p and used with the exact density p(r), the agreement 
seen in the above figures would be completely spoiled. 
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FIG. 15: Test of the TF functional relation ([T02|l for = 
100068 particles in the three-dimensional spherical billiard 
(lines as in Fig. 1141 units ft^/2m = R = 1; both densities 
divided by A''^''''). Note that in the vertical direction of the 
figure, only a very small excerpt around the bulk value is 
displayed. 



VI. SUMMARY AND CONCLUDING 
REMARKS 

We have presented a semiclassical theory, initiated in 
for the oscillating parts of the spatial densities in 
terms of closed non-periodic orbits (NPOs), while the 
smooth part of the densities are given by the (extended) 
Thomas-Fermi (TF) theory. Our equations (gOl) - (O 
are the analogues of the semiclassical trace formula ([55)) 
for the density of states in terms of periodic orbits. 

For spherical systems, two kinds of oscillations in the 
spatial densities can be separated, as is implied in Eqs. 
([^ - ([^ : regular, short-ranged ones (denoted by the 
symbol Sr) that we can attribute to the librating NPOs 
in the radial direction, and irregular, long-ranged ones 
(denoted by 5i„) that are due to non-radial NPOs and 
therefore only exist in _D > 1 dimensions. The simple 
nature of the radial NPOs leads immediately to a number 
of relations between the regular parts of the oscillations, 
such as Eqs. ([50]) . ([55)) . or the universal form ([^ for 
(5rp(r) valid near r = 0. It also explains that the kinetic- 
energy density ^(r) defined in © has no rapid regular 
oscillations, as implied in ()26p . but is smooth for all one- 
dimensional systems, as well as for isotropic harmonic 
oscillators [iCf and linear potentials [1^ in arbitrary D 
dimensions, since these contain no non-radial NPOs. 

In spherical systems, the semiclassical expansion in 
terms of NPOs is expected to work best for filled "main 
shells" where the total energy has a pronounced local 
minimum. This is also discussed in Ref. [2l| on the two- 
dimensional circular billiard, for which a complete clas- 
sification of all NPOs (in addition to the periodic orbits) 
has been made and the semiclassical theory for the spatial 
density oscillations has been studied analytically. The 
semiclassical approximation for the density oscillations 
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is, indeed, found there to work best for the closed-shell 
systems with filled main shells. But even for "mid-shells" 
systems with half-filled main shells and for most interme- 
diate systems, the agreement of the semiclassical densi- 
ties with the quantum-mechanical ones has turned out in 
[21I I to be very satisfactory. 

Based on the semiclassical theory, we were able to gen- 
eralize the "local virial theorem" (LVT) given in (j49|l and 
(jlOip . which had earlier been derived from exact results 
for isotropic harmonic oscillators [l^l, to arbitrary local 
potentials V{r). We emphasize that the LVT is valid 
(semiclassically) also for an interacting N -fermion sys- 
tem bound by the self-consistent Kohn-Sham potential 
obtained within the framework of DFT and might be 
verified experimentally in finite fermionic systems. 
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APPENDIX: INCLUSION OF FINITE 
TEMPERATURE IN THE SEMICLASSICAL 
THEORY 

In this Appendix we give a short sketch of how to in- 
clude finite-temperatures in the semiclassical formalism. 
Extensions of semiclassical trace formulae to finite tem- 
peratures have been used already long ago in the context 
of nuclear physics [60t and more recently in mesoscopic 
physics [IJ . We shall present here a derivation by means 
of a suitable folding function, which has proved useful 
also in the corresponding microscopic theory [6l| . 

For a grand-canonical ensemble of fermions embedded 
in a heat bath with fixed temperature, the variational 
energy is the so-called grand potential 51 defined by [1^ 



Fermi occupation numbers 



n= (H) -TS - X(N) 



(A.l) 



where H and N are the Hamilton and particle number 
operators, respectively, T is the temperature in energy 
units (i.e., we put the Boltzmann constant ks equal to 
unity), S is the entropy, and A the chemical potential. 
Note that both energy and particle number are conserved 
only on the average. For non-interacting particles, we can 
write the Helmholtz free energy F as 



F={H) - TS 



(A.2) 



where £"„ is the energy spectrum of H and Un are the 



1 



1 + eXp jT 

and the entropy is given by 



(A.3) 



5 = -2 ^ K logl^n + (1 - Vn) l0g(l - Vn)] ■ (A.4) 
n 

The chemical potential A is determined by fixing the av- 
erage particle number 



TV = (TV) =2^, 



(A.5) 



It can be shown [6l| that the above quantities TV, F 
and S may be expressed in terms of a convoluted finite- 
temperature level density gxiE) as 



F = 2 / EgT{E)AE. 

J —00 



(A.6) 



The function griE) is defined by a convolution of the 
"cold" (T = 0) density of states ([13| 



9t{E) = r g{E') ME - E') dE' ^ ^ ^E - K) , 

J -co „ 

(A.7) 



whereby the folding function friE) is given as 
ME) ^ 



ATcos\i{E/2T) 



(A., 



Note that all sums in (|A.2[) - (|A.8[) run over the complete 
(infinite) spectrum of the Hamiltonian H. It is now easily 
seen that 

N = 2 I gT{E)dE. (A.9) 



To show that the integral (|A.6p gives, indeed, the correct 
free energy (jA.2p including the "heat energy" —TS needs 
some algebraic manipulations. From F, the entropy S 
can always be gained by the canonical relation 



S 



dF 



(A.IO) 



The same convolution can now be applied also to the 
semiclassical trace formula ([M]) for the oscillating part of 
the density of states which we re-write as 



Sg{E) ~ RcY,Apo(E)e-k 



Spo(-E)— iCTpo 



(A.ll) 



PO 



The oscillating part SgxiE) of the finite-temperature 
level density is obtained by the convolution of (jA.lip 
with the function friE) as in (jA.7|) . In the spirit of 
the stationary-phase approximation, we take the slowly 
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varying amplitude ^po (E) outside of the integration and 
approximate the action in the phase by 

5po {E') ~ 5po (E) + {E' -E)Tpo{E), (A. 12) 

so that the resuh becomes a modified trace formula 

6gT{E) ~ Re^^Po(£;)/T(rpo(i?))e*^-°(^)-'^^°, 

(A.13) 



PO 



where 



Tpo{E)^Tpo{E)/h 



(A. 14) 



and the temperature modulation factor is given by 
the Fourier transform of the convolution function fx'. 



The Fourier transform of the function (jA.Sp is known, so 
that 



Mr) 



ttTT 



sinh(7rTT) 



(A.16) 



The "hot" trace formula (jA.lSp with the modulation fac- 
tor (|A.16p has been obtained in [i^, [l^l • 

For the spatial densities, we can proceed exactly in the 
same way. For the particle density, e.g., the microscopic 
expression ([2]) is replaced by 



PT(r) = 2^|0„(r)|2z.„, 



(A.17) 



where the sum again runs over the complete spectrum. 
Starting from the semiclassical expression (|l(7)l for Sp{r) 
at T = 0. we rewrite it as 



Spo{X, r) ~ Re ^ ^npo(A, r) e'* 



(A.18) 



NPO 



where the amplitude ^npo collects all the prefactors of 
the phase in ([40]) . The finite-T expression is given by the 
convolution integral 



(5/9t(A, r) 



dpo{X-E,r)fTiE)dE. (A.19) 



Expanding the phase under the integral as above, we 
arrive at 



(A. 15) 6pT{X,r)^RcJ2-^ 

NPo(A, r) /t(%po(A, r)) 



(A,r) 



NPO 

(A.20) 

where Tnpo = Tnpo(A, r)//i is the period of the NPO in 
units of h. The corresponding expressions for the other 
spatial densities are obvious. 

For the smooth parts of the densities, we recall that 
the (E)TF theory at T > is well known and refer to 
chapter 4.4.3 of [IHl for the main results and relevant 
literature. 

Other types of correlations can be included in the semi- 
classical theory in the same way, as soon as a suitable 
folding fcori {E) function corresponding to friE) in I A.8I) 
and its Fourier transform are known (see. e.g., Ref. 63|). 
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TABLE I: Contributions of different types of non-periodic orbits to the spatial densities in a rectangle billiard with sides Qo; 
and Qy. The first row gives the position of the images of P with {k^, ky) G . The second row gives the length of the orbit 
and the third row the angle 9 between the initial and final momentum. The fourth and fifth rows give the contributions to 5p 
and 5t\, respectively. 

Type of orbits a b c 

Image points of P(a::,y) {2ka:Qx + x,2kyQy — y) {2kxQx — x,2kyQy + y) {2kxQx — x,2kyQy — y) 

Orbit length L{kxQx, kyQy — y) L{kxQx — x,kyQy) L{kxQx — x,kyQy — y) 

9^ = -2 arctan ( ''"^JfJ ^ 9^ = 2 arctan ^ klcit-x 1 ^ c = tt 
Contribution to 5p Spi, = f{kxQx,kyQy — y,l) Spb = f{kxQx — x,kyQy,l) 5pc = f{kxQx — x, kyQy — y,0) 

Contribution to Sti 5ria — Xcos{9^)5pa Stu, — Acos(Sb)'5pb Stic — \cos{9c)Spc 



